search.arima <- function(x, d=NA, D=NA, max.p=5, max.q=5,
    max.P=2, max.Q=2, max.order=5, stationary=FALSE, ic=c("aic","aicc","bic"),
    trace=FALSE,approximation=FALSE,xreg=NULL,offset=offset,allowdrift=TRUE,
    parallel=FALSE, num.cores=NULL)
{
    #dataname <- substitute(x)
    ic <- match.arg(ic)
    m <- frequency(x)
         
    oldwarn <- options()$warn
    options(warn=-1)
    on.exit(options(warn=oldwarn))

    if(allowdrift)
        maxK <- (d+D <= 1)
    else
        maxK <- ((d+D) == 0)

    # Choose model orders
    #Serial - technically could be combined with the code below
    if (parallel==FALSE)
    {
        best.ic <- 1e20
        for(i in 0:max.p)
        {
            for(j in 0:max.q)
            {
                for(I in 0:max.P)
                {
                    for(J in 0:max.Q)
                    {
                        if(i+j+I+J <= max.order)
                        {
                            for(K in 0:maxK)
                            {
                                fit <- myarima(x,order=c(i,d,j),seasonal=c(I,D,J),constant=(K==1),trace=trace,ic=ic,approximation=approximation,offset=offset,xreg=xreg)
                                if(fit$ic < best.ic)
                                {
                                    best.ic <- fit$ic
                                    bestfit <- fit
                                    constant <- (K==1)
                                }
                            }
                        }
                    }
                }
            }
        }
    } else
################################################################################
    # Parallel
    if (parallel==TRUE){

		to.check <- WhichModels(max.p, max.q, max.P, max.Q, maxK)
		
            par.all.arima <- function(l){
                .tmp <- UndoWhichModels(l)
                i <- .tmp[1]; j <- .tmp[2]; I <- .tmp[3]; J <- .tmp[4]; K <- .tmp[5]==1

                if (i+j+I+J <= max.order){
                    fit <- myarima(x,order=c(i,d,j),seasonal=c(I,D,J),constant=(K==1),trace=trace,ic=ic,approximation=approximation,offset=offset,xreg=xreg)
                }
                if (exists("fit")){
                    return(cbind(fit, K))
                } else return(NULL)
            }
			
			if(is.null(num.cores)) {
				num.cores <- detectCores()
			}
			
            # clusterApplyLB() for Windows, mclapply() for POSIX
            if (Sys.info()[1] == "Windows"){
                cl <- makeCluster(num.cores)
                all.models <- clusterApplyLB(cl=cl, x=to.check, fun=par.all.arima)
                stopCluster(cl=cl)
            } else all.models <- mclapply(X=to.check, FUN=par.all.arima, mc.cores=num.cores)

            # Removing null elements
            all.models <- all.models[!sapply(all.models, is.null)]

            # Choosing best model
            best.ic <- 1e20
            for (i in 1:length(all.models)){
                if(!is.null(all.models[[i]][, 1]$ic) && all.models[[i]][, 1]$ic < best.ic){
                    bestfit <- all.models[[i]][, 1]
                    best.ic <- bestfit$ic
                    constant <- unlist(all.models[[i]][1, 2])
                }
            }
            class(bestfit) <- "Arima"
        }
################################################################################
    if(exists("bestfit"))
    {
        # Refit using ML if approximation used for IC
        if(approximation)
        {
            #constant <- length(bestfit$coef) - ncol(xreg) > sum(bestfit$arma[1:4])
            newbestfit <- myarima(x,order=bestfit$arma[c(1,6,2)],
                seasonal=bestfit$arma[c(3,7,4)],constant=constant,ic,trace=FALSE,approximation=FALSE,xreg=xreg)
            if(newbestfit$ic > 1e19)
            {
                options(warn=oldwarn)
                warning("Unable to fit final model using maximum likelihood. AIC value approximated")
            }
            else
                bestfit <- newbestfit            
        }
    }
    else
        stop("No ARIMA model able to be estimated")

    bestfit$x <- x
    bestfit$series <- deparse(substitute(x))
    bestfit$ic <- NULL
    bestfit$call <- match.call()
    #bestfit$call$data <- dataname
#    bestfit$xreg <- xreg

    if(trace)
        cat("\n\n")

    return(bestfit)
}


ndiffs <- function(x,alpha=0.05,test=c("kpss","adf","pp"))
{
    test <- match.arg(test)
    require(tseries)
    x <- c(na.omit(c(x)))
    d <- 0
    oldwarn <- options(warn=-1)
    if(test=="kpss")
        dodiff <- kpss.test(x)$p.value < alpha
    else if(test=="adf")
        dodiff <- adf.test(x)$p.value > alpha
    else if(test=="pp")
        dodiff <- pp.test(x)$p.value > alpha
  else
    stop("This shouldn't happen")
    if(is.na(dodiff))
    {
        options(warn=oldwarn$warn)
        return(d)
    }
    while(dodiff & d<2)
    {
        x <- diff(x)
        d <- d+1
        if(test=="kpss")
            dodiff <- kpss.test(x)$p.value < alpha
    else if(test=="adf")
      dodiff <- adf.test(x)$p.value > alpha
    else if(test=="pp")
      dodiff <- pp.test(x)$p.value > alpha
    else
      stop("This shouldn't happen")
        if(is.na(dodiff))
            return(d-1)
    }
    options(warn=oldwarn$warn)
    return(d)
}

# Set up seasonal dummies using Fourier series
SeasDummy <- function(x)
{
    n <- length(x)
    m <- frequency(x)
    if(m==1)
        stop("Non-seasonal data")
    tt <- 1:n
    fmat <- matrix(NA,nrow=n,ncol=2*m)
    for(i in 1:m)
    {
        fmat[,2*i] <- sin(2*pi*i*tt/m)
        fmat[,2*(i-1)+1] <- cos(2*pi*i*tt/m)
    }
    return(fmat[,1:(m-1)])
}

# CANOVA-HANSEN TEST
# Largely based on uroot package code for CH.test()
SD.test <- function (wts, s=frequency(wts))
{
    if(any(is.na(wts)))
        stop("Series contains missing values. Please choose order of seasonal differencing manually.")
    if(s==1)
        stop("Not seasonal data")
    t0 <- start(wts)
    N <- length(wts)
    if(N <= s)
        stop("Insufficient data")
    frec <- rep(1, as.integer((s+1)/2))
    ltrunc <- round(s * (N/100)^0.25)
    R1 <- as.matrix(SeasDummy(wts))
    lmch <- lm(wts ~ R1,na.action=na.exclude)   # run the regression : y(i)=mu+f(i)'gamma(i)+e(i)
    Fhat <- Fhataux <- matrix(nrow=N, ncol=s-1)
    for (i in 1:(s-1))
        Fhataux[, i] <- R1[,i] * residuals(lmch)
    for (i in 1:N)
    {
        for (n in 1:(s - 1))
            Fhat[i, n] <- sum(Fhataux[1:i, n])
    }
    wnw <- 1 - seq(1, ltrunc, 1)/(ltrunc + 1)
    Ne <- nrow(Fhataux)
    Omnw <- 0
    for (k in 1:ltrunc)
        Omnw <- Omnw + (t(Fhataux)[, (k + 1):Ne] %*% Fhataux[1:(Ne - k), ]) * wnw[k]
    Omfhat <- (crossprod(Fhataux) + Omnw + t(Omnw))/Ne
    sq <- seq(1, s-1, 2)
    frecob <- rep(0,s - 1)
    for (i in 1:length(frec))
    {
       if (frec[i] == 1 && i == as.integer(s/2))
           frecob[sq[i]] <- 1
       if (frec[i] == 1 && i < as.integer(s/2))
           frecob[sq[i]] <- frecob[sq[i] + 1] <- 1
    }
    a <- length(which(frecob == 1))
    A <- matrix(0, nrow=s - 1, ncol=a)
    j <- 1
    for (i in 1:(s - 1)) if (frecob[i] == 1)
    {
        A[i, j] <- 1
        ifelse(frecob[i] == 1, j <- j + 1, j <- j)
    }
    stL <- (1/N^2) * sum(diag(solve(t(A) %*% Omfhat %*% A, tol=1e-25) %*% t(A) %*% t(Fhat) %*% Fhat %*% A))
    return(stL)
}


forecast.Arima <- function (object, h=ifelse(object$arma[5] > 1, 2 * object$arma[5], 10),
    level=c(80, 95), fan=FALSE, xreg=NULL, lambda=object$lambda, ...)
{
#    use.constant <- is.element("constant",names(object$coef))
    use.drift <- is.element("drift", names(object$coef))
    if (is.element("x", names(object)))
        x <- object$x
    else
        x <- object$x <- eval.parent(parse(text=object$series))
    usexreg <- (!is.null(xreg) | use.drift | is.element("xreg",names(object)))# | use.constant)
#    if(use.constant)
#        xreg <- as.matrix(rep(1,h))
    if(!is.null(xreg))
    {
        xreg <- as.matrix(xreg)
        h <- nrow(xreg)
    }
    if (use.drift)
    {
        n <- length(x)
        if(!is.null(xreg))
            xreg <- cbind((1:h)+n,xreg)
        else
            xreg <- as.matrix((1:h)+n)
    }
    if(usexreg)
    {
        if(is.null(xreg))
          stop("No regressors provided")
        if(!is.null(object$xreg))
            object$call$xreg <- object$xreg
        else # object from arima() rather than Arima()
        {
            xr <- object$call$xreg
            object$call$xreg <- if (!is.null(xr))         
                                    eval.parent(xr)
                                else NULL
        }
        if(ncol(xreg) != ncol(object$call$xreg))
          stop("Number of regressors does not match fitted model")
        pred <- predict(object, n.ahead=h, newxreg=xreg)
    }
    else
        pred <- predict(object, n.ahead=h)

    # Fix time series characteristics if there are missing values at end of series.
    tspx <- tsp(x)
    nx <- max(which(!is.na(x)))
    if(nx != length(x))
    {
        tspx[2] <- time(x)[nx]
        start.f <- tspx[2]+1/tspx[3]
        pred$pred <- ts(pred$pred,f=tspx[3],s=start.f)
        pred$se <- ts(pred$se,f=tspx[3],s=start.f)
    }
        
    if(fan)
        level <- seq(51,99,by=3)
    else
    {
        if(min(level) > 0 & max(level) < 1)
            level <- 100*level
        else if(min(level) < 0 | max(level) > 99.99)
            stop("Confidence limit out of range")
    }

    nint <- length(level)
    lower <- matrix(NA, ncol=nint, nrow=length(pred$pred))
    upper <- lower
    for (i in 1:nint)
    {
        qq <- qnorm(0.5 * (1 + level[i]/100))
        lower[, i] <- pred$pred - qq * pred$se
        upper[, i] <- pred$pred + qq * pred$se
    }
    colnames(lower)=colnames(upper)=paste(level, "%", sep="")
    method <- arima.string(object)
    fits <- fitted(object)
    if(!is.null(lambda))
    {
      pred$pred <- InvBoxCox(pred$pred,lambda)
      lower <- InvBoxCox(lower,lambda)
      upper <- InvBoxCox(upper,lambda)
  }
    return(structure(list(method=method, model=object, level=level,
        mean=pred$pred, lower=lower, upper=upper, x=x,
        xname=deparse(substitute(x)), fitted=fits, residuals=residuals(object)),
        class="forecast"))
}


forecast.ar <- function(object,h=10,level=c(80,95),fan=FALSE, lambda=NULL, ...)
{
    pred <- predict(object,n.ahead=h)
    if(fan)
        level <- seq(51,99,by=3)
    else
    {
        if(min(level) > 0 & max(level) < 1)
            level <- 100*level
        else if(min(level) < 0 | max(level) > 99.99)
            stop("Confidence limit out of range")
    }
    nint <- length(level)
    lower <- matrix(NA,ncol=nint,nrow=length(pred$pred))
    upper <- lower
    for(i in 1:nint)
    {
        qq <- qnorm(0.5*(1+level[i]/100))
        lower[,i] <- pred$pred - qq*pred$se
        upper[,i] <- pred$pred + qq*pred$se
    }
    colnames(lower)=colnames(upper)=paste(level,"%",sep="")
    method <- paste("AR(",object$order,")",sep="")
    x <- as.ts(eval.parent(parse(text=object$series)))
    f=frequency(x)
    res <- ts(object$resid[-(1:object$order)],start=tsp(x)[1]+object$order/f,f=f)
    fits <- x-res
    
    if(!is.null(lambda))
    {
      pred$pred <- InvBoxCox(pred$pred,lambda)
      lower <- InvBoxCox(lower,lambda)
      upper <- InvBoxCox(upper,lambda)
      fits <- InvBoxCox(fits,lambda)
      x <- InvBoxCox(x,lambda)
    }

    return(structure(list(method=method,model=object,level=level,mean=pred$pred,lower=lower,upper=upper,
        x=x, xname=deparse(substitute(x)), fitted=fits,residuals=res)
        ,class="forecast"))
}

# Extract errors from ARIMA model (as distinct from residuals)
arima.errors <- function(z)
{
  if(!is.list(z))
    stop("z must be a list")
  if(is.element("x",names(z)))
    x <- z$x
  else
  {
    series.name <- z$series
    if(is.null(series.name))
      stop("missing component series in argument z\n")
    x <- eval.parent(parse(text=series.name))
  }
  if(!is.element("xreg",names(z)))
  {
    if(!is.element("xreg",names(z$coef)))
      return(x)
    else
      xreg <- eval.parent(z$coef$xreg)
  }
  else
    xreg <- z$xreg
  norder <- sum(z$arma[1:4])
  if(is.element("intercept",names(z$coef)))
    xreg <- cbind(rep(1,length(x)),xreg)
  return(ts(x - xreg %*% as.matrix(z$coef[(norder+1):length(z$coef)]),f=frequency(x),s=start(x)))
}

# Return one-step fits
fitted.Arima <- function(object,...)
{
  if(is.element("x",names(object)))
    x <- object$x
  else
  x <- eval.parent(parse(text=object$series))
  if(is.null(object$lambda))
  return(x - object$residuals)
  else
    return(InvBoxCox(BoxCox(x,object$lambda) - object$residuals, object$lambda))
}

# Calls arima from stats package and adds data to the returned object
# Also allows refitting to new data
# and drift terms to be included.
Arima <- function(x, order=c(0, 0, 0),
      seasonal=list(order=c(0, 0, 0), period=NA),
      xreg=NULL, include.mean=TRUE, include.drift=FALSE, include.constant, lambda=model$lambda,
    transform.pars=TRUE,
      fixed=NULL, init=NULL, method=c("CSS-ML", "ML", "CSS"),
      n.cond, optim.control=list(), kappa=1e6, model=NULL)
{
    # Remove outliers near ends
    #j <- time(x)
    #x <- na.contiguous(x)
    #if(length(j) != length(x))
    #    warning("Missing values encountered. Using longest contiguous portion of time series")

  series <- deparse(substitute(x))
  
  origx <- x
  if(!is.null(lambda))
    x <- BoxCox(x,lambda)

  if (!is.null(xreg))
  {
    nmxreg <- deparse(substitute(xreg))
    xreg <- as.matrix(xreg)
    if (is.null(colnames(xreg))) 
        colnames(xreg) <- if (ncol(xreg) == 1) nmxreg else paste(nmxreg, 1:ncol(xreg), sep="")
  }
  
  if(!missing(include.constant))
  {
    if(include.constant)
    {
      include.mean <- TRUE
      if((order[2] + seasonal$order[2]) == 1)
        include.drift <- TRUE
    }
  }
    
  if(!is.null(model))
  {
    tmp <- arima2(x,model,xreg=xreg)
    xreg <- tmp$xreg
  }
  else
  {
    if(include.drift)
    {
      drift <- 1:length(x)
      xreg <- cbind(drift=drift,xreg)
    }
    if(is.null(xreg))
      tmp <- stats:::arima(x=x,order=order,seasonal=seasonal,include.mean=include.mean,
          transform.pars=transform.pars,fixed=fixed,init=init,method=method,n.cond=n.cond,optim.control=optim.control,kappa=kappa)
    else
      tmp <- stats:::arima(x=x,order=order,seasonal=seasonal,xreg=xreg,include.mean=include.mean,
             transform.pars=transform.pars,fixed=fixed,init=init,method=method,n.cond=n.cond,optim.control=optim.control,kappa=kappa)
  }
  tmp$series <- series
  tmp$xreg <- xreg
  tmp$call <- match.call()
  tmp$lambda <- lambda
  tmp$x <- origx
  
  return(tmp)
}

# Refits the model to new data x
arima2 <- function (x, model, xreg)
{
    use.drift <- is.element("drift",names(model$coef))
    use.intercept <- is.element("intercept",names(model$coef))
    use.xreg <- is.element("xreg",names(model$call))
    if(use.drift)
    {
      driftmod <- lm(model$xreg[,"drift"] ~ I(time(model$x)))
      newxreg <- driftmod$coeff[1] + driftmod$coeff[2]*time(x)
      if(!is.null(xreg))
        xreg[,"drift"] <- newxreg
      else
        xreg <- as.matrix(data.frame(drift=newxreg))
      use.xreg <- TRUE
    }
   
    if(model$arma[5]>1 & sum(abs(model$arma[c(3,4,7)]))>0) # Seasonal model
    {
        if(use.xreg)
            refit <- Arima(x,order=model$arma[c(1,6,2)],seasonal=list(order=model$arma[c(3,7,4)],period=model$arma[5]),
                fixed=model$coef,include.mean=use.intercept,xreg=xreg)
        else
            refit <- Arima(x,order=model$arma[c(1,6,2)],seasonal=list(order=model$arma[c(3,7,4)],period=model$arma[5]),
                fixed=model$coef,include.mean=use.intercept)
    }
    else if(length(model$coef)>0) # Nonseasonal model with some parameters
    {
        if(use.xreg)
           refit <- Arima(x,order=model$arma[c(1,6,2)],fixed=model$coef,xreg=xreg,include.mean=use.intercept)
        else
            refit <- Arima(x,order=model$arma[c(1,6,2)],fixed=model$coef,include.mean=use.intercept)
    }
    else # No parameters
            refit <- Arima(x,order=model$arma[c(1,6,2)],include.mean=FALSE)

    refit$var.coef <- matrix(0,length(refit$coef),length(refit$coef))
    if(use.xreg) # Why is this needed?
      refit$xreg <- xreg
    return(refit)
}

# Modified version of function in stats package

print.Arima <- function (x, digits=max(3, getOption("digits") - 3), se=TRUE,
    ...)
{
#    if (!is.element("x", names(x)))
#        x$x <- eval.parent(parse(text=x$series))

    cat("Series:",x$series,"\n")
    cat(arima.string(x),"\n")
  if(!is.null(x$lambda))
    cat("Box Cox transformation: lambda=",x$lambda,"\n")
    #cat("\nCall:", deparse(x$call, width.cutoff=75), "\n", sep=" ")
#    if(!is.null(x$xreg))
#    {
#        cat("\nRegression variables fitted:\n")
#        xreg <- as.matrix(x$xreg)
#        for(i in 1:3)
#            cat("  ",xreg[i,],"\n")
#        cat("   . . .\n")
#        for(i in 1:3)
#            cat("  ",xreg[nrow(xreg)-3+i,],"\n")
#    }
    if (length(x$coef) > 0) {
        cat("\nCoefficients:\n")
        coef <- round(x$coef, digits=digits)
        if (se && nrow(x$var.coef)) {
            ses <- rep(0, length(coef))
            ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits=digits)
            coef <- matrix(coef, 1, dimnames=list(NULL, names(coef)))
            coef <- rbind(coef, s.e.=ses)
        }
        print.default(coef, print.gap=2)
    }
    cm <- x$call$method
    if (is.null(cm) || cm != "CSS")
    {
        cat("\nsigma^2 estimated as ", format(x$sigma2, digits=digits),
            ":  log likelihood=", format(round(x$loglik, 2)),"\n",sep="")
        npar <- length(x$coef) + 1
        nstar <- length(x$residuals) - x$arma[6] - x$arma[7]*x$arma[5]
        bic <- x$aic + npar*(log(nstar) - 2)
        aicc <- x$aic + 2*npar*(nstar/(nstar-npar-1) - 1)
        cat("AIC=", format(round(x$aic, 2)), sep="")
        cat("   AICc=", format(round(aicc, 2)), sep="")
        cat("   BIC=", format(round(bic, 2)), "\n",sep="")
    }
    else cat("\nsigma^2 estimated as ", format(x$sigma2, digits=digits),
        ":  part log likelihood=", format(round(x$loglik, 2)),
        "\n", sep="")
    invisible(x)
}

# Modified version of function in stats package

predict.Arima <- function(object, n.ahead=1, newxreg=NULL, se.fit=TRUE, ...) 
{
    myNCOL <- function(x) if (is.null(x)) 
        0
    else NCOL(x)
    rsd <- object$residuals
    ## LINES ADDED
    if(!is.null(object$xreg))
        object$call$xreg <- object$xreg
    ## END ADDITION
    xr <- object$call$xreg
    xreg <- if (!is.null(xr)) 
        eval.parent(xr)
    else NULL
    ncxreg <- myNCOL(xreg)
    if (myNCOL(newxreg) != ncxreg) 
        stop("'xreg' and 'newxreg' have different numbers of columns: ", ncxreg, " != ", myNCOL(newxreg))
    class(xreg) <- NULL
    xtsp <- tsp(rsd)
    n <- length(rsd)
    arma <- object$arma
    coefs <- object$coef
    narma <- sum(arma[1:4])
    if (length(coefs) > narma) {
        if (names(coefs)[narma + 1] == "intercept") {
            xreg <- cbind(intercept=rep(1, n), xreg)
            newxreg <- cbind(intercept=rep(1, n.ahead), newxreg)
            ncxreg <- ncxreg + 1
        }
        xm <- if (narma == 0) 
            drop(as.matrix(newxreg) %*% coefs)
        else drop(as.matrix(newxreg) %*% coefs[-(1:narma)])
    }
    else xm <- 0
    if (arma[2] > 0) {
        ma <- coefs[arma[1] + 1:arma[2]]
        if (any(Mod(polyroot(c(1, ma))) < 1)) 
            warning("MA part of model is not invertible")
    }
    if (arma[4] > 0) {
        ma <- coefs[sum(arma[1:3]) + 1:arma[4]]
        if (any(Mod(polyroot(c(1, ma))) < 1)) 
            warning("seasonal MA part of model is not invertible")
    }
    z <- KalmanForecast(n.ahead, object$model)
    pred <- ts(z[[1]] + xm, start=xtsp[2] + deltat(rsd), frequency=xtsp[3])
    if (se.fit) {
        se <- ts(sqrt(z[[2]] * object$sigma2), start=xtsp[2] + 
            deltat(rsd), frequency=xtsp[3])
        return(list(pred=pred, se=se))
    }
    else return(pred)
}
